!> author: 左志华
!> date: 2022-07-02
!>
!> Internal forces <br>
!> 内力：粘性力、压力
!> @todo 找到解决初始时刻边界压力和排斥力扰动问题
module sph_internal_force

    use sph_kinds, only: rk
    use sph_particle, only: particle_t
    use twin_sph_list_type, only: twin_list_t
    use sph_eos, only: p_ideal_gas, p_water
    use sph_env, only: nums

    procedure(internal_force_fcn), pointer :: internal_force => internal_force_with_visc
                                                                        !! Internal force function pointer <br>
                                                                        !! 内力函数指针
    logical :: viscosity    !! Viscosity on/off <br>
                            !! 是否考虑粘性力

    abstract interface

        !> Internal force function <br>
        !> 内力函数
        subroutine internal_force_fcn(particle, pairs, acc, n)
            import :: particle_t, twin_list_t, rk
            type(particle_t), intent(inout) :: particle  !! 区域
            type(twin_list_t), intent(in) :: pairs  !! 粒子对
            real(rk), intent(inout) :: acc(:, :)  !! 加速度
            integer, intent(in) :: n  !! 区域内总粒子数
        end subroutine internal_force_fcn

    end interface

contains

    !> Internal force function without viscosity <br>
    !> 内力求解
    pure subroutine internal_force_without_visc(particle, pairs, acc, n)
        type(particle_t), intent(inout) :: particle  !! 区域
        type(twin_list_t), intent(in) :: pairs  !! 粒子对
        real(rk), intent(inout) :: acc(:, :)  !! 加速度
        integer, intent(in) :: n  !! 区域内总粒子数
                                  !! @note 未实际使用
        integer :: i
        real(rk) :: h(2)

        ! ----------------------------- 根据 EoS 计算压强信息 -------------------------------- !
        do concurrent(i=1:nums(4))

            select case (particle%itype(i))
            case (-1, 1); call p_ideal_gas(particle%rho(i), particle%u(i), particle%p(i), particle%c(i))  ! 理想气体
            case (-2, 2); call p_water(particle%rho(i), particle%p(i), particle%c(i))             ! 淡水
            end select

        end do

        ! ------------------------------------ 计算压力 ------------------------------------ !
        do concurrent(i=1:pairs%log)
            associate (ip => pairs%items(1, i), jp => pairs%items(2, i))
                if (ip > nums(1) .or. jp > nums(1)) cycle            ! 双虚粒子不参与内力计算

                ! 用到了虚粒子的压强对实粒子的影响，防止边界处实粒子有较大的、指向边界的初始压强合力
                h(:) = -(particle%p(ip) + particle%p(jp))*pairs%dwdx(:, i)/(particle%rho(ip)*particle%rho(jp))
                acc(:, ip) = acc(:, ip) + h(:)*particle%mass(jp)
                acc(:, jp) = acc(:, jp) - h(:)*particle%mass(ip)
            end associate
        end do

    end subroutine internal_force_without_visc

    !> Internal force function with viscosity <br>
    !> 内力求解（含粘性力）
    !> @todo 粘性熵、温度变化
    pure subroutine internal_force_with_visc(particle, pairs, acc, n)
        type(particle_t), intent(inout) :: particle  !! 区域
        type(twin_list_t), intent(in) :: pairs  !! 粒子对
        real(rk), intent(inout) :: acc(:, :)  !! 更新加速度
        integer, intent(in) :: n  !! 实粒子总数量
                                  !! @note 用以声明本地数组长度
        integer :: i
        real(rk), parameter :: r23 = 2.0_rk/3.0_rk  !! 2/3
        real(rk) :: h(2), dvel(2), hxx, hxy, hyy
        real(rk) :: txx(n), txy(n), tyy(n), eta(n)

        txx = 0.0_rk; txy = 0.0_rk; tyy = 0.0_rk    ! 初始化本地数组

        ! ----------------------------- 根据 EoS 计算压强等信息 -------------------------------- !
        do concurrent(i=1:nums(4))

            select case (particle%itype(i))
            case (-1, 1)    ! 理想气体
                call p_ideal_gas(particle%rho(i), particle%u(i), particle%p(i), particle%c(i))
                eta(i) = 0.0_rk
            case (-2, 2)    ! 淡水
                call p_water(particle%rho(i), particle%p(i), particle%c(i))
                eta(i) = 1.0e-3_rk  ! 运动学粘性系数
            end select

        end do

        do concurrent(i=1:pairs%log)
            associate (ip => pairs%items(1, i), jp => pairs%items(2, i))
                if (ip > nums(1) .or. jp > nums(1)) cycle         ! 双虚粒子不参与内力计算

                ! ----------------------------- 计算剪切张量 ---------------------------------- !
                dvel(:) = particle%vel(:, jp) - particle%vel(:, ip)
                hxx = r23*(2*dvel(1)*pairs%dwdx(1, i) - dvel(2)*pairs%dwdx(2, i))
                hxy = r23*(dvel(1)*pairs%dwdx(2, i) + dvel(2)*pairs%dwdx(1, i))
                hyy = r23*(2*dvel(2)*pairs%dwdx(2, i) - dvel(1)*pairs%dwdx(1, i))

                txx(ip) = txx(ip) + particle%mass(jp)*hxx/particle%rho(jp)
                txx(jp) = txx(jp) + particle%mass(ip)*hxx/particle%rho(ip)
                txy(ip) = txy(ip) + particle%mass(jp)*hxy/particle%rho(jp)
                txy(jp) = txy(jp) + particle%mass(ip)*hxy/particle%rho(ip)
                tyy(ip) = tyy(ip) + particle%mass(jp)*hyy/particle%rho(jp)
                tyy(jp) = tyy(jp) + particle%mass(ip)*hyy/particle%rho(ip)

                ! ----------------------------- 计算压力、剪切力 ------------------------------ !
                ! 用到了虚粒子的压强对实粒子的影响，防止边界处实粒子有较大的、指向边界的初始压强合力
                ! 采用 \( (p_i + p_j)/(rho_i*rho_j) \) 方式
                h(:) = -(particle%p(ip) + particle%p(jp))*pairs%dwdx(:, i)
                ! he = sum((particle%vel(:, jp) - particle%vel(:, ip))*h(:))
                h(1) = h(1) + &
                   &(eta(ip)*txx(ip) + eta(jp)*txx(jp))*pairs%dwdx(1, i) + &
                    (eta(ip)*txy(ip) + eta(jp)*txy(jp))*pairs%dwdx(2, i)
                h(2) = h(2) + &
                   &(eta(ip)*txy(ip) + eta(jp)*txy(jp))*pairs%dwdx(1, i) + &
                    (eta(ip)*tyy(ip) + eta(jp)*tyy(jp))*pairs%dwdx(2, i)

                associate (rho_ij => particle%rho(ip)*particle%rho(jp))
                    acc(:, ip) = acc(:, ip) + h(:)*particle%mass(jp)/rho_ij
                    acc(:, jp) = acc(:, jp) - h(:)*particle%mass(ip)/rho_ij
                end associate

            end associate
        end do

    end subroutine internal_force_with_visc

end module sph_internal_force
